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Essential glossary 



• analytical: precise, exact. Analytical methods are based on analytical in- 
version formulas. 

• attenuation: physical phenomenon of the deviation of an electromagnetic 
ray from its natural path, due to the contact of a body. 

• CT: X-ray Computed Tomography, a medical imaging procedure that uti- 
lizes computer-processed X-rays to produce 'slices' of specific areas of the 
body. 

• dummy: (also called phantom, we use that name to avoid confusions) 
physical device that simulates a certain organ and used to make experiments 
on the tomography machine, in vitro. 

• electromagnetic waves: is an oscillation which propagates in the space, 
transporting energy but not matter. We distinguish different kind of waves 
by their wavelength, among which: 

— X-rays: wavelength in the range of 0.01 to 10 nanometers. 

— Gamma-rays: wavelength shorter than 1 picometer. 

• ill-conditioned problem: given f(x) = y, y € Y, we have to find x G X, 
with X and Y vector spaces. This problem is i.e. if an arbitrarily small 
perturbation on y can lead to an arbitrarily large error on x. 

• in silico: experiment made on mathematical simulation of the data (in this 
work we use phantoms) and not using real data. 

• in vitro: experiments made on real data coming from a biological/chemical 
phenomenon reproduced in a test tube (in this work we use a dummies), 
and not using an organism. 

• kernel function: a function K such that K(x, y) = <ft(x — y) with (ft a RBF. 

• nuclear medicine: a branch of medicine involving the application of ra- 
dioactive substances in the diagnosis and treatment of disease. 



• 



over(/under)-determined: a problem with more(/less) equations than 
unknowns. Given f(x) = y, y £ Y, we have to find x £ X, with X and 
Y vector spaces. This problem is o.d. if dim(Y) > dim(X) and u.d. if 
dim(Y) < dim(X). 

phantom: an image created to test the efficiency of the algorithms in silico, 
without need of data given by a physical machine. The most used phantom 
in this work is the Shepp-Logan phantom, which is similar to a brain, see 
Figure 3.2. 

Randon transform: is the integral transform consisting of the integral of a 
function over straight lines, proposed by the Austrian mathematician Johann 
Radon in 1917. Is simulates mathematically the effect of attenuation of 
some electromagnetic waves passing through different materials at several 
angles, and then it is useful for the CT reconstruction. 

RBF: Radial Basis Function, functions ^ : H C I s — )■ M. which depends 
only on the radial component r = \\x\\. 

sampling: times or points at which the data is collected from a continuous 
phenomenon. 

scattering: erroneous acquisition of some rays made by the tomography 
machine due to attenuation. 

sinogram: anlxp matrix of data collected by the tomography machine in 
I angular scans. In the case of in silico experiments we call s. the discrete 
Radon transform of the phantom. 

SPECT: Single-Photon Emission Computed Tomography, is a nuclear medicine 
imaging technique using radioactive tracers that emit Gamma-rays in de- 
caying. 
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Introduction 



)j£"?W™^(\ he problem of image reconstruction from projections dates back 
c^\ B/^ ^° ^ e 1970s w ith the first X-ray machines, which give out as 
■^2yHfe<^ data the result of some circular acquisitions. An algorithm 
f^^zM!^^ provides to extract the relevant section of the body, as if seen 
from above. We discuss about these algorithms and their efficiency both 
in the case of the transmission tomography (such as the well-known CAT, 
which emits and detects X-rays), and of the emission tomography (for in- 
stance PET, which detects pairs of positrons and SPECT, which captures 
gamma-rays) used in nuclear medicine 1 . In the first case a series of rays 
pass through the section before arriving at a detector in order to obtain a 
section of the organs, while in the second case the rays are coming from in- 
side the patient's body, in which they were previously injected tracers (not to 
be confused with the radiocontrast agents) in order to obtain cross-sectional 
areas in which is concentrated and the tracer, e.g. an area affected by tumor 
or a working area of a certain organ. 

In particular, we will discuss the hybrid SPECT/CT case, in which the 
SPECT image, obtained with a radioactive tracer with low emissions, is re- 
lated to a low-resolution CT scan on the same section, in order to understand 
the relation between the area highlighted by the marker and the neighboring 
organs. 



1 nuclear medicine uses radioactive tracers, i.e. chemical elements with low radioactivity 
which are injected into the patient and which the human body uses as if they were normal 
atoms from physiological reactions. E.g. the ls FdG, 18-Flu-deoxy-glucose used in PET, 
which the body absorbs as if it were glucose, of which the tumor cells are "greedy". From 
the annihilation or from the decay of these tracers is possible to recover their position. 
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Chapter 1 

Some medical and technical issues 



Mais vous, mon cher Pangloss, dit Candide, comment se 
peut-il que je vous revoie? II est vrai, dit Pangloss, que 
vous m'avez vu pendre; je devais naturellement etre brule 
mais vous vous souvenez qu'il plut a verse lorsqu'on allait 
me cuire: l'orage fut si violent qu'on desespera d'allumcr 
le feu: je fus pendu, parcequ'on ne put mieux faire: un 
chirurgien acheta mon corps, m'emporta chez lui, et me 
dissequa. II me fit d'abord une incision cruciale depuis le 
nombril jusqu'a la clavicule. On ne pouvait pas avoir ete 
plus mal pendu que je l'avais ete. [..] Enfin je respirais 
encore: l'incision cruciale me fit jeter un si grand cri, que 
mon chirurgien tomba a la renverse. 

Voltaire, Candide 



n this chapter we will discuss some medical, physical and techni- 
cal issues useful to better understanding the image reconstruc- 
tion problem and to decide what assumptions we need in order 
to model the behavior of the machines with sufficient accuracy 
and thus to obtaining a good reconstruction. 



1.1 Electromagnetic waves 

An electromagnetic wave is an oscillation which propagates in the space, 
transporting energy but not matter. Its intuitive form is a cosine function 
with amplitude A and frequency v. Its speed in the vacuum is c = 3 • 10 8 m/s, 
the speed of light, so if A is the wavelength, then \v = c. The energy carried 




by the wave is E — hv (where h = 6 • 10 u Js is the Plank constant), hence 
the energy is inversely proportional to the wavelength: 

Some examples of electromagnetic waves are the light of the sun, infrared, 
microwaves, radio waves, up to X-rays and gamma rays. 

It is known as X-rays that portion of the electromagnetic spectrum with 
a wavelength in the range of 10~ 3 up to 10 1 nm discovered in 1895 by the 
German physicist Wilhelm Conrad Rontgen, who earned for this discover the 
first Nobel price in 1901. X-rays are produced from the transition of elec- 
trons, and they can penetration through many objects, included the human 
body. For this reason X-rays are widely used in medicine to discover irregu- 
larities into the patient's body. 

Similar to X-rays, but shorter in wavelength (then they transport more 
energy for the (1.1)), the gamma rays (sometimes written as 7-rays) are 
produced by a radioactive decaying 1 . As the X-rays, gamma rays can be 
dangerous and cause burns, cancer, and genetic alteration, so doctors must 
pay attention not to do useless or too hazardous exams, by limiting the 
energy of the waves used and the time of exposition. 



1.2 Attenuation 

Every electromagnetic waves passing thought a material mean is subject to 
an attenuation due to diffraction; in every point the ray has a probability 
p to deviate from its natural path 2 . The linear attenuation coefficient /i is 
proportional to the probability of deviation while passing in 1cm of matter; 
/i depends on the material, on its density and on the energy carried by the 
wave. The linear attenuation coefficient is proportional to the density p, so 
the mass attenuation coefficient fi m = -, only depends on the energy and the 
material. 

Table 1.1 shows the mass attenuation coefficients of different materials (wa- 
ter, sodium iodide and lead) and energies. 



Hhis is the fundamental difference between X and gamma rays: the first ones come 
from outside the nucleus of the atoms, the second ones from inside. 

2 see [2] Ch 6D, p. 80 for a complete and detailed survey about photon attenuation and 
scattering. 



Photon energy 




Mass attenuation coefficient fi m 




(eV) 




(cm 2 /g) 




H 2 


Nal 


Pb 


80 


0.179 


3.00 


2.07 


100 


0.168 


1.64 


5.23 


150 


0.149 


0.590 


1.89 


200 


0.136 


0.314 


0.945 


300 


0.118 


0.158 


0.383 


400 


0.106 


0.112 


0.220 


500 


0.097 


0.092 


0.154 


600 


0.089 


0.080 


0.120 


800 


0.079 


0.066 


0.085 


1000 


0.070 


0.058 


0.069 



Table 1.1: Mass attenuation coefficient for several values of the photon energy and 
for different materials. 

It is evident from this table that the mass attenuation coefficient is de- 
creasing with the growth of the energy or with the decreasing of the density, 
but it is neither inversely proportional to the energy or proportional to the 
density. Another popular way to represent this attenuation is the Hounsfield 
scale 

jl fJ'water 
(J"water 

which compares the attenuation coefficient of the medium and that of the 

water. 

Computerized tomography, in its classical form, calculates the attenuation 

coefficient of the different areas inside the body from the difference between 

entering and outgoing X-rays, in order to recognize the different organs and 

find possible malfunctions or malformations. 



1.3 Scattering 

During the reconstruction we have to face different kind of issues related to 
the attenuation. The most known, called scattering, is the effect of a mul- 
tiple attenuation of the same ray. In fact assume that a ray is diffracted in 
the point (xi,yi) and deviates from its natural path; assume also that, as 
in Figure 1.1, it deviates again in the point (£2,2/2) m the direction of the 
detector (perpendicular in the figure). This means that the machine will see 
many rays "out of place", but with low energies. 



The most common way to reduce the scattering effect 3 is a Hamming win- 
dowing of the incoming data around the known peak of energy of the used 
tracer or ray, as shown in Figure 1.2. 



outgoing ray 




Counts 



Figure 1.1: Scattering effect. 



20% window 




120 140 160 



Detected 
Source 



Energy (KeV) 



Figure 1.2: Scattering correction example for 99m Tc, who has a single peak of 
energy of 140KeV. 



3 See [2] Ch IOC p. 161 for further explanations. 



1.4 Partial volume effect 



Another typical issue in tomography is the correct reconstruction of the 
higher- valued areas. In fact in such area the reconstruction will result smoother 
and lower-valuer than the original ones. Conversely the lower-valued areas 
next to it will have a growth of value near the higher area. This is known as 
partial volume effect, and it is similar to the light bulb effect: it is hard 
to see the exact position of of the bulb since it is on because of the light rays 
next to it. As we will see later the partial volume effect is not related to the 
reconstruction algorithm we use, but only depends from resolution. 




Original shape 

Reconstruction 



Figure 1.3: Partial volume effect example on the characteristic function of an 
interval. 



1.5 Structure of a Gamma camera 




Figure 1.4: External view of a Gamma-camera. 

A machine for axial tomography externally has a toroidal structure. Inside of 
it is placed the patient through a sliding bed. The position of the couch ob- 
viously determines the height of the body of the patient at which we want to 
obtain the section. Some algorithms give a three-dimensional reconstruction 
from the data taken at different heights; other a two-dimensional reconstruc- 
tion for each set of data relating to the same height. Inside the machine 
there are the emitters (only in the case of a transmission tomography) and 
the detectors. 
Now, let us describe the structure of the detectors of the Gamma-camera 



6 



used in SPECT; the structure of other types of detectors are similar. 
The gamma rays leaving the highlighted area are selected by the collima- 
tors. The collimators are devices with circular, square or hexagonal section 
long L, and placed parallel to each other at a constant distance and nor- 
mal to the circumference of capture. The hole has length D and the septal 
thickness (the space between two holes) is t. Their aim is to filter incoming 
rays passing only those nearly parallel (this is named parallel-hole tomogra- 
phy). Other dispositions of the collimators give different data acquisition 
(divergent/convergent-hole or pin-hole tomography). 
The ray that passes out the collimators arrives in a scintillator crystal that 



energy discrimination 
and positioning electronics 



s( 



E 



t D 




photomultiplier 
tubes 



^J scintillator 
crystal 



collimators 



highlited area 



Figure 1.5: Structure of a Gamma-camera used in SPECT. 



converts the beam into a spark of low-density This spark gets in the pho- 
tomultiplier tubes (or PM tubes), which converts the sparkle into variable 
beams of electrons with an energy of O(10 5 ) times greater than that of the 
sparkle. The machine will successively treat these beams as electrical pulses, 
convert to numerical data, and detect the location of the pulses; this last 
process is made by an algorithm which tries to indicate from which hole be- 
tween the collimators had passed the original gamma ray. 
Furthermore we call c the distance between the collimators and crystal, x 



the object-collimator distance and L e ff — L the effective length of the 

collimator with \x attenuation coefficient of the material of the collimators 4 . 
From this knowledge on the functioning of the machine we can derive the 
sensitivity and the resolution of the instrument 5 . Once estimated the resolu- 
tion and sensitivity, we also get the a 'priori error, which we will use during 
the experiments. 

1.6 Resolution and sensitivity 

1.6.1 Septal thickness 

Does the collimator really stop all the non-parallel 7-rays? There is not any 
material able to stop all the rays, and, by consequence we have to set a level 
of septal penetration which can be considered negligible (e.g. a tolerance 
of 5 percent) and set the thickness as the corresponding value. Let w be 
the shortest path length for 7-rays to travel from one hole to the next as in 
Figure 1.6. Then, since t <^ I, we find the proportion t + 2D : I — t : w, 
which leads to the equality w = ^ D . 
Since the septal penetration has to be less than 5 percent, then for Beer- 




crystal 



collimators 



7-ray 



Figure 1.6: Scheme for the minimal septal thickness which gives the desired toler- 
ance; note that t + 2D : I = t : w. 



4 usually lead, /1 = 2.49mm *. 

5 For further details on the following equations see [10] Ch 5, p. 87. 



Lambert's law 6 e ttw < 0.05; we know that e 3 = 0.05, then w > -, thus 

^ 6D 

t> 



/il — 3 

This means that if we want to reduce the septal thickness the collimators must 
be made of material with high density and atomic number. Lead is the most 
used material because is cheap and easy to find; other materials, as tantalum, 
tungsten, and even gold have been employed in experiments 7 . From the 
formula above we also deduce that if the energies emitted by the source is 
bigger, then the attenuation coefficient will decrease, and the minimal septal 
thickness must be larger. 

1.6.2 Geometric resolution 

What is the smallest distance between two objects such that the machine can 
recognize them as two? To answer this question let us consider a point source 
of gamma rays as in Figure 1.5, and the two red rays in the figure (the parallel 
one and the last entering one). If R is the distance of the two rays at the 
PM tubes, then, by the similitude of triangles L e ff : D = (L e ff + x + c) : R, 
hence 

r = d(i + x + " 



L eff 

since c<x the formula becomes 

R = D + x 



L 



eff 



The resolution of the machine depends on the collimator resolution and on 
the intrinsic resolution (the resolution of the crystal and the electronics) 



So we have to find the intrinsic resolution to find the resolution of the system. 

1.6.3 Experimental resolution 

Another, equivalent way to find the resolution of the machine is to scan a 
little sphere filled of a certain quantity of tracer (a point source). The Full 
Width at Half Maximum (briefly FWHM) is the distance between the points 



6 it will be introduced in the next chapter. 

7 See [17] Ch 13, p. 332 for a more detailed explanation. 
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Figure 1.7: The FWHM of reconstructed data from a sphere filled with tracer. 
Note that is similar to a Gaussian function. 



of level |, where h is the maximum of the function. E.g. if the function is 



the Gaussian f(x) = L- exp — 2o .2 it is easy to see that its FWHM is 

2-^/2 log(2) cr^ 2.355 a. 

The FWHM of the reconstructed data 8 offers a realistic, direct measure of 

the resolution of the machine. In fact two point sources distant less than 

FWHM are recognized as a unique, bigger source. On the other hand if their 

distance is greater than FWHM we can see the contribution of two different 

sources see Figure 1.8. 

The resolution has also a role in partial volume effect: if the resolution is 

smaller, partial volume effect are less intense, because the sampling with a 

machine of a certain FWHM it is equivalent to the windowing with Gaussian 

function of a = FWHM : the windowed function will be smoother if a is 



bigger. 



FWHM . 
2-^/21(^(2)' 



1.6.4 Sensitivity 

The sensitivity of the detectors is given by the ratio between gamma ray 
passed through the collimators and all the gamma rays emitted by the source, 
i.e. is proportional to the solid angle subtended by irradiated detector area 

8 as we will see in the last chapter, the reconstruction method and the filtering of the 
image can change our estimate of the FWHM. In particular a low-pass filter increases the 
partial volume effect and the resolution rises by consequence. 
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counts 




FWHM FWHM 



pixels 



Figure 1.8: The FWHM is the limit for two objects to be recognizable as two; 
the sum of the two signals (dashed) at a distance R can be recognized as the 
contribution of two different objects. 

multiplied by the square of the fraction of a collimator covered by the hole: 

<kR 2 D 2 

b oc 



,2 x 2 



S = K 2 



D 2 



where K is a constant given by the shape of the collimator section (K = 0.28 
for square holes, 0.26 for hexagonal, 0.24 for round). Since D ^> t the formula 
becomes 



S 



KD 



L 



eff 



All these formulas will be used during the experiments. 

1.7 Tracers used in SPECT and standard ex- 
ams 

The tracers used in SPECT are depending on the body part under examina- 
tion and of the type of analysis to be performed. 

Table 1.2 shows the main types of studies that use SPECT, and reports the 
respective tracer with its half-life, and some data about the projections and 
the size of the image that we want to get. 
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Type of Study 


Substance 


Emission 
Energy 


Half-Life 


#Angles 


Size Image 


Projection 
Time 


Bone Scan 


99mrn 


140KeV 


6h 


120 


128x128 


30s 


Myocardial 
perfusion scan 


99mrn 


140KeV 


6h 


60 


64x64 


25s 


Brain scan 


99mrp 


140KeV 


6h 


64 


128x128 


30s 


Neuroendocrine 
or neurological 
tumor scan 


123j /131j 


159KeV 


13h/8days 


60 


64x64 


30s 



Table 1.2: Main types of studies that use SPECT tomography. 

From this table we see that the size of the image to be reconstructed and 
the number of angles of sampling varies, but the typical values are 128x128 
for the image size (or 64x64) and 120 (or 60) for the number of projections. 
Furthermore the time of acquisition for each angle is about 30s, hence the 
maximum exposure time is 60 minutes (that is also the limit of average 
patient's tolerance to stay inside the machine), which is far less than the 
half-life of tracers. For this we can assume that during the acquisition energy 
emitted by the tracer remains constant, i.e. does not depend on an acquisition 
time when it is conducted. 

Let us take a look to the tracers. The 99m Tc technetium-99m is a nuclear 
isomer of the "Tc ("m" stands for "metastable"); its half life is 6 hours, so 
the 93% of the atoms decay in 24 hours; moreover its radioactive emission is 
made up almost exclusively of gamma radiation with a single peak of energy 
(140KeV), sufficient for a good penetration through the body structures, but 
not dangerous for the patient. The 99m Tc has the ability to bind chemically 
in many biologically active molecules, e.g. when is combined with tin it binds 
blood cells and can therefore be used to discover eventual disorders of the 
circulatory system. 

123 J and 131 J are radioisotopes of Iodine; they are used in radiotherapy, to 
cure hyperthyroidism and some types of thyroid cancer that absorb iodine. 
They emit gamma and beta rays so they can be seen in diagnostic scans after 
its use as therapy. 

1.8 Some assumptions 

In order to model correctly the tomography machine we will make some 
assumptions: 
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1. rays (X, gamma, etc..) All have no width, i.e. their length is much 
greater than their width, so we can consider them as straight lines. 

2. incoming rays and detected rays are parallel to the line of angle 9 
(parallel-beam); other cases must be specified. 

3. the machine runs / scans at the angles (6?i, ...,#;); the detectors move 
on a circle (the angles are typically a uniform partition of [0, 2tt] or 
[0,tt] with a step of 3°or 6°). 

4. if the total time of acquisition is relatively short is indifferent to run 
an acquisition before or after another one 9 . 

5. the reconstructed image has dimension N x N (typically iV = 128 or 
N = 64). 

6. for each scan j G 1, . . . , / the machine gets p linear samples (sij, . . . , s p j) 
(typically p = N, or p = 185 or p = 95). So the sinogram is an / x p 
matrix of data collected in / scans each of p elements. 

7. in the case of transmission tomography the ray which passes through 
the body is subject to an unknown attenuation f(x,y), of which we 
know the sum along some lines at different angles. 

8. in the case of emission tomography the rays start from the inside of 
the body and is subject to attenuation a(x, y) supposed to be note. 
Our goal is locating the radioactive concentration f(x,y), given the 
attenuated sum along some lines at different angles. 



''see previous section. 

13 



Chapter 2 

Mathematical modeling: the 
Radon transform and its variants 



Nella vita del signor Palomar c'e stata un'epoca in cui la 
sua regola era questa: primo, costruire nella sua mente 
un modello, il piu perfetto, logico, geometrico possibile; 
secondo, verificare se il modello s'adatta ai casi pratici 
osservabili nell'esperienza; terzo, apportare le correzioni 
necessarie perche modello e realta coincidano. 

I. Calvino, Palomar 

n this chapter different mathematical models of the tomography 
problem are given as transforms of functions, each coming from 
different assumptions on the projection method and on the kind 
of tomography performed. 



2.1 The Beer-Lambert's law and the Radon trans- 
form 

The Radon transform, introduced in 1917 by the Austrian mathematician 
Johann Karl August Radon, can be physically interpreted as a ray passing 
through different materials, each one of them with its own attenuation coeffi- 
cient. We know the intensity in which the rays entered in and released by the 
body, so we can find the attenuation coefficients. In fact, Beer-Lambert's law 
says that if I(x) is the intensity of the ray, A(x) the attenuation coefficient 
at the point x, satisfies 

AI = -A(x)I(x) Ax. 
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Equivalently, by integrating their respective differential terms, 

C^-Ct-MS) (2 - l) 

where the second term is known. Equation (2.1) has an equivalent on the 
2-dimensional plane. In fact, given the change of coordinates due to rotation 

x \ _ / cos 9 — sin 9 \ ( t 
y J y sin 9 cos 9 I V s 

we get the formulation of the Radon transform as an integral on a line 

Kf(t,6)= [ f= [ f(x)5(t-x-9)dx 

Je (t:6) Jr 2 

where x — (x,y), 9 — (cos 9, sin 9), or equivalently: 

Kf(t,6)= J f(tcos9-ssm9,tsm9 + scos9)ds= [ f \t9 + s6 L ) ds (2.2) 

Jr Jr 

where 9 X = (—sin 9, cos 9). 

2.2 Alternative Transforms 

2.2.1 Backprojection 

Now we may ask for every point (x, y) what is the average of the rays passing 
through that point. To answer this question we introduce the formal adjoint 
operator to TZ, the so-called backprojection operator 

TZ*g(x ) y) = —-— \ g(x cos 9+y sin 9, 9) d9 = — — / g(s : 9)S(s— x-9) dsd9 

\b I Js 1 P I JrxS 1 

(2.3) 
or equivalently 

n*g(x,y) = j^ [ g(x-9,9)d9 
P I Js 1 

where S* 1 is the unit hemisphere represented by the parameter 9, i.e. S 1 = 
[0, 7r] or S 1 = [0, 2n] according to the angular sampling, so IS* 1 ) = 7r or 2n. 
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Figure 2.1: Radon transform of an object with constant attenuation coefficient. 



2.2.2 Divergent beam transform 

In the Radon transform formulation we assumed that the collimators receive 
only the parallel rays. Rays are not parallel, but start from a single trans- 
mitter and are disposed radially. If now we let pass through the collimators 
all the rays, we have to consider the fan beam transform, also called di- 
vergent beam transform 

/■+00 r+oo 

Vh(x,9)= h(x + tcos6,y + tsm9)dt = h(x + W)dt (2.4) 

Jo Jo 

where the parameter t this time indicates the distance traveled by the ray 
from the emitter. 



2.2.3 Attenuated Radon transform 

In the case of an emission tomography, as PET or SPECT, modeling the 
tomography becomes more difficult: in fact each single beam emitted by 
the tracer in the highlighted area must pass through objects of different 
consistency, and then will be subject to attenuation. From Beer-Lambert's 
law we obtain that 



I(xi) = J(x )exp 



■fi 



A(x)dx 
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We assume then to know the attenuation coefficients 1 a(x,y), we know that 
a ray passing through (x, y) changes direction with a certain probability 
directly proportional to a(x,y). Then, if f(x,y) is the concentration of the 
tracer, then the radiation detected at an angle 9 to a first approximation is 
given by 

K a f(t, 9) = [ e~ a ^f(x) = n(e~ a f)(t, 9). (2.5) 

which is simply the Radon transform of a (weighted) function e~ a f. 

We can obtain a more realistic version of (2.5) if we consider that a ray 
subjected to diffraction expands radially from the point of diffraction 

TZJ(t,9)= [ e- Va W + Vf(x) (2.6) 

This operator is called attenuated Radon transform (briefly AtRT), where 
T> is the divergent beam transform. 

Note that this new transform is not equivalent to the Radon transform of a 
modified of weighted function. 

In the special case of PET tomography, where the positrons come out 
two by two in opposite directions, the attenuation part in (2.6) goes outside 
the integral 2 and then to recover the section we only have to apply an at- 
tenuation correction to the image obtained by the inversion of the classic 
Radon transform to obtain a realistic reconstruction. If the scan is done on 
a relatively homogeneous areas such as the brain, we can also consider the 
attenuation as a constant 3 . 



Hhe attenuation coefficients were previously unknown and called f(x,y), now we call 
them a(x,y) and we suppose to know them, while the unknown function is now f(x,y), 
the tracer's concentration 

2 see [13] Ch 1.2, p.4, and[18] Ch 6, p.91 for more details. 

3 more details on the dedicated chapter. 
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Chapter 3 

Analytical properties and 
inversion of the Radon transform 



[...] as if all this had been done and said so many times 
before it made you feel it was recorded, they all in here 
existed basically as Fourier Transforms of postures and 
little routines, locked down and stored and call-uppable 
for rebroadcast at specified times. 

D.F. Wallace, Infinite Jest 

et us show some properties of the Radon transform and its 
related transforms; some of this properties will be useful in 
^!> the following chapters, some other will lead us to an analytical 
resolution formula to the following problem 

Given g find / such that IZf = g (3-1) 

Let us begin with some definitions: 

Definition 1. A function f : M 2 — > K. is called (Radon-) transformable 
iflZf(t, 9) exists and is finite V(t, 8) G R x S 1 . 

We will not go deeper in discussing the domain of the Radon transform, 
but it's good to know from [9] and [13] that the Radon transform can be seen 
both as a map from the Schwartz space S in a subset of S, or from L 2 to 
itself or from C^° to itself 1 . Let us define in general 

Definition 2. A function f is transformable (as regards one of the opera- 
tors defined in previous chapter) if its transform exists and is finite (almost) 
everywhere. 

^^most of theorems shown below require / £ S or / £ C^°. 
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3.1 Basic properties 

Lemma 3.1.1. Let f be a Radon-transformable function. Hence 

Rf(t,e)=Kf(-t,ir + $) 

Proof. Trivial by substitution in (2.2) of page 16. □ 

We can therefore restrict the range of angular projections in 9 G [0,7r], 
as the projections from the angles from n to 2n do not provide new data; on 
the contrary in the case of attenuated transform the sampling must be made 
of angles from to 2-k. 



Lemma 3.1.2. The operators 1Z, 1Z* , T>, lZ a are linear on their respective 
spaces of transformable functions. 

Proof. This is also trivial by the linearity of the integral operator. E.g. 
K(af + 0g)(t, 9) = [ {af{x) + (3g(x))6(t - x ■ 9) dx = 



n 



a I f(x)S(t -x-9) dx + /3 / g(x)8(t - x ■ 9) dx = a Kf(t, 9) +/3 TZg{t, 

'R 2 



Proposition 3.1.3. For every appropriated functions f and g, then 

K(f* 2 g) = Kf*Kg 
n*(g*Kf)=K*g*1Zf 

where "*" is meant to be the convolution made only on the linear variable t 
and "*2 " the 2-dimensional convolution. 

Proof, see [13] Ch 2.1, p. 13-14 for the proof. □ 

Next Lemmas will be useful in interpolation. 

Lemma 3.1.4. If f is transformable and radial, i.e f(x) = (p(\\x\\), then IZf 
only depends on the linear parameter t 



Proof. 



Kf(t,< 



'y 



f(t9 + s9 ± )ds= / <p(\\t9 + s9 ± \\)ds 



cos U — sin 6 
sin# cos^ 



p(Vt 2 + s 2 ) ds 



because the matrix is orthogonal. 



D 
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Theorem 3.1.5. (Shift property) Let c be a constant vector and g(x) 
f(x — c). Then 

TZg{t,e) = nf(t-c-e,e) 



Proof. See [16] Ch 2.7, p.31. D 

This last Theorem guarantees the uniqueness of the reconstruction prob- 
lem (3.1) in the continuous case. 

Theorem 3.1.6. (Uniqueness) Let f E S and K C IR 2 be a convex, com- 
pact set called the hole. IflZf = in every plane {x ■ 9 = s} fl K = 0. 
Then f = outside K . 

Proof. Proof can be found in [13] Ch 2.3, p.30-32. □ 



3.2 Sampling and resolution 

In order to solve the problem of reconstruction (3.1) in the discrete case we 
have to make some further assumptions on the function /. We shall also 
make some assumptions on the sampling to avoid aliasing problems. 
Let us begin with a definition and a standard Theorem in applied mathe- 
matics 

Definition 3. A Fourier-transformable function f is called b-band-limited 

if its Fourier transform J-f(v) is negligible for \v\ > b, i.e. 

\Ff(a)\ da = 

\u\>b 

Theorem 3.2.1. (Shannon-Nyquist) Let f be a b-band-limited function 
and let h < ir jb Then f is uniquely determined by the values f(hk), fcGZ" 
and in L 2 (W a ) 

,7T 



f(x) = y^ f(hk)sinc(—(x — hk)) 

k 

L2 = h^(j2\f(hk)\ 2 ) 1/2 



-h 

k 

moreover we have 



Proof, the proof is shown in [13] Ch 3.1, p. 56. □ 
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This theorem gives a condition to recover efficiently a 6-band-limited func- 
tion: the sampling step h must be smaller than ir/b (this is called Nyquist 
condition). 

Now we want to find a condition on the angular step for recover in a reli- 
able way the function /. Let us proceed with some further definitions and 
properties. 

Definition 4. A spherical harmonic function of degree m is the restric- 
tion to S 1 of harmonic, homogeneous polynomial function of degree m on 

n 2 

Proposition 3.2.2. If H m is the space of the spherical harmonic functions 
of degree < m inM. 2 then a basis for this space is 

{PO,P(l,l),P(l,-l)}l=l,...,m 

with 

Po = 1 

P(,,1)(0)=COS(W) 
?(,,_!)($) =Sill(J0) 

where, as usual 9 = (cos 9, sin 9) 

We can easily verify that P(i,±i) are harmonic and homogeneous polyno- 
mials in 9; in fact, by de Moivre's formula 

cos(nx) + i sin(nx) = (cos x + i sin x) n 

we can easily find the explicit functions by taking 

P(i t ±i) = 9lc/3m(cos 6 + % sin 9) 1 

developed with the Newton's binomial formula. Then these functions are all 
homogeneous polynomials in (cos 9, sin 9) = (x,y). On the other side these 
functions are harmonic too, since their laplacians are 

A (p(i,+i)) + iA (p(i,-i)) = A (p(i,+i) + m-i)) = A ((^ + w) n ) = 

= n(n - l)(x + iy) n ~ 2 + {i) 2 n(n - l)(x + iy) n ~ 2 = 
Now we define a subset of this space which will be useful for next theorems 

Definition 5. The space H' m is the subset of the functions in H m which are 
even function for m even and odd functions for m odd. 
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We observe immediately that dim(H' m ) — to + 1 because the basis of this 
space is {p(j,(_i)i)}j=o...m- 

Definition 6. A C S 1 is m-resolving if the only function on H' m vanishing 
in A is the function p = 0. 

Lemma 3.2.3. If A is m-resolving then \A\ > dim(H' m ) — to + 1. 

Theorem 3.2.4. Let A be a m-resolving set and f E C£°(B(0, 1)). IflZf 
vanishes in A, then for < d < 1 

ll/IU-< * -M f,$m) 

27rr]{v,m) 

where rj is a function such that 

< 77(1?, to) < C(i?) exp(-A(i?)TO) 
and 6 > B(d), B, C, A positive, and 

e o (f,0m) = ±- [ \Tg(a,e)\da 

Zn J\a\>i)m 

Proof, see [13] Ch 3.3, p. 68 for the proof and for further details. □ 

This last theorem gives an answer to the question of resolution. In fact if 
f is 6-band-limited then e is negligible and / can be recovered reliably from 
a m-resolving set, where m > b, i.e. where the angular sampling points are 
at least / = m + 1 = \b\ + 2. 

Definition 7. Let f E C^°(S(0, 1)) be b-band limited and such that IZf 
vanishes on a l-resolving set. A sampling of IZf of the form {(£«,#•/)} of 
linear step h and I angular data, is said that IZf is reliably sampled if 
I — 1 > b and h < ir/b. 

In conclusion if we take / angular samples of the projection data each one 
made by p linear equispaced acquisitions of step h = t—t, then the maximum 
frequency of the data will be 

b = min{/ - 2, ^-p} (3.2) 

and the resolution will be 

-[^sampling ■ 

with attention to report the radius of sampling to r instead of 1, i.e. we have 
to multiply R samv iing by r. 
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3.3 Inversion formulas for the Radon transform 

Proposition 3.3.1. The backprojection operator is the Hermitian added to 
the transform 1Z, but not its inverse: in fact 



K*Kf(x) = 2 



m 



R2 f - y\ 



-Av 



7^*2 f 

\x\ 



Proof, see [13] Ch 2.1, p. 15 for the proof. 



□ 



this means that what we get from applying the backprojection operator 
to projection data is a blurred image. From this property we can derive a 
first numerical method of inversion which simply uses the 2D FFT and IFFT 
algorithms for a fast deconvolution. The computational complexity of such 
algorithm is known to be 0(N 2 log N) where N is the number of pixel by 
size in the reconstructed image. 



Theorem 3.3.2. (Central Slice Theorem) Let f be a Fourier and Radon 
transformable function. Then 

T 2 f(T cos 6, T sin 6) = 7 (Rf) (T, 9) 

Where we mean for Ti the 2-dimensional Fourier transform and J 7 the 1- 
dimensional Fourier transform applied only to the variable t. 

Proof, proof can be found in [5] Ch 6, p. 47; a more general result can be 
found in [13] Ch 2.1, p.10. □ 



n 



m 




Figure 3.1: Scheme for Central Slice Theorem: the 2-dimensional Fourier transform 
of a function in polar coordinates is equal to the 1-dimensional Fourier transform 
made on the linear parameter of its Radon transform. 
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Theorem 3.3.3. (Inverse Radon transform) The original function can 
be derived analytically using the filtered backprojection formula 

f = \K* [^(MJW))] (3-3) 

where we mean that the direct and inverse Fourier transform is applied 
only to the variable t. 

Proof, see [5] Ch 6, p. 49 or [18] Ch 7, p. 46 for the proof. □ 

In what follows we will refer this formula as FBP. 

Given the inaccurate numerical results of this formula, we use 

w{v) = v(u)\u\ 

instead of \i/\, with v a low-pass filter 2 . The reconstruction formula becomes 

/ = \lV [F-\w{v)HKf))\ (3.4) 

Sometimes we know w(t) = T~ 1 {w{i/)) in order to apply the convolution 
directly to the IZf, without going through the Fourier transforms, i.e. 

/ - l -K* [w{t) * Kf] 

An example of filter is the function w(u) = XulliIH w ith L fixed. 

3.4 Discretization of the FBP formula and er- 
ror bound 

We can get another approach proposed by [13] in Ch 5.1 for the resolution of 
the problem (3.1). In fact thanks to Theorem (3.1.3), if W b = lZ*(w b ), then 

W b *f = K*(w b *Kf) 

if we assume that ui} } is such that W b approximates the 5— distribution we get 

f = n*(w b *nf) (3.5) 

which is equivalent to (3.4). As said in the previous section we will choose 
W b as a low-pass filter with b its cut-off frequency. 

2 for a list of the low-pass filter used see the chapter about regularization. 
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Original phantom f(x,y) 



sinogram: Radon trans orm Rf(s ,th eta) 





Filtered backprojection 



Unflltered backprojection 





Figure 3.2: Original Shepp-Logan phantom, sinogram of its Radon transform, and 
reconstruction by filtered and unflltered back projection formulas. 
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Now we consider the discrete version of (3.1). Let g = IZf be sampled at 

(ti, 6j) with i = —q, . . . , q and j = 1, . . . , / and let h = 1/q be the linear 

step, i.e. ti = hi; in this notation p — 2q + 1. If we consider the discrete 

convolution 

q 

a*b(t) = h^ a(t - ti)b(ti) 

i=-q 

and the discrete backprojection given from a quadrature rule 

1 P 

71 3=1 

with cx(pj) positive weights, now the formula (3.5) reads as 

Ifbp = n;(w b 'inf) (3.6) 

The following theorem gives an error estimate for this last formula. 

Theorem 3.4.1. Let f G C£°(B(0, 1)), assume that the quadrature rule used 
before is exact in H' 2m , and let IZf be reliably sampled, i. e., for some d G (0, 1) 



then 
and 



b < <dm b < n/h 



f FBP = W b * f + ei + e 2 



N<^(/,») 



where 



JcJ,6)\da 



|e 2 |<77(t?,m)||/|Uo 

£*(/,&) = 1^1 sup / \Tg{c 

ees 1 J\a\>b 

and f] is any function such that 

< 77(1?, to) < (7(i?) exp(-A(i?)m) 
and b > B($), B, C, A positive. 
Proof, see [13] Ch 5.1, p. 104 for the proof. □ 
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Let us try to explain the meaning of this theorem. If / is 6-band-limited, 
the term e\ is negligible, so the first error, related to the inaccuracy of the 
discretization of the convolution product, is due to aliasing problems. On 
the other side the second error term e 2 depends on the chosen discretization 
and the discretization of the backprojection operator. 

Moreover from the proof 3 of Theorem (3.4.1) we can find a weaker result 
which will be more useful for our purposes 

|e 2 | < 2|S' 1 | \\w h *g-(w h *g) m \\ L oo 

where g = IZf and (wb*g) m is the truncation to order m = b of the spherical 
harmonics series. Since iuj, is given exactly we can assume 

w b * g - (w b * g) m = w b * (g -g) 

with g the truncation to order b of the spherical harmonics series for g, i.e. 
g — g is the error on the approximation of g — TZf. 



Lemma 3.4.2. For every suitable functions a and b on a suitable set Q C W, 

s GN, we get 

||o * &||l°° < HflllLiH^IlL 00 

Proof. 

||a*6|| L oo = sup | / a(x-y)b(y) dy\ < sup / \a(x - y)\\b(y)\ dy < 

< sup \b(x)\ sup / \a(x — y)\ dy — sup \b(x)\ \ \a(y)\ dy — ||o||x,i||6||i,oo 

X X Jq X Jq 

a 

Remark: this is a particular case of Young's inequality for convolutions. 

We are now able to prove the following Theorem 

Theorem 3.4.3. Let f E C£°(B(0, 1)) a b-band-limited function, and let 
g = IZf be reliably sampled, i.e. 

b < m and b < n/h 

let f be the FBP reconstruction, then 

\f(x) - f(x)\ < 2|S' 1 | ||wb||Li||fi'-^||L- + |e 3 | 
3 this proof is available in [13], as said previously. 
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Vx 6 .8(0,1), i.e. 

11/ - /||l°°(r 2 ) < 2\S I ||w6||x,i(R)||5 - glU^Rxs 1 ) + |e 3 j 

w;#/j e3 the quadrature error which depends on the interpolation chosen to 
approximate the backprojection integral. 

We will see some application of this Theorem in the Chapter of the ex- 
periments. 
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Chapter 4 

Algebraic Reconstruction 
Techniques 



Insomnia il secondo tO in cui stanno la freccia FO e un 
po' piu in la il leone LO e qui il me stesso QO e uno strato 
spaziotemporale che resta fermo e identico per sempre, e 
accanto ad esso si dispone tl con la freccia Fl e il leone 
LI e il me stesso Ql che hanno leggermente cambiato le 
loro posizioni, e li affiancato c'e t2 che contiene F2 e L2 e 
Q2 e cosi via. In uno di questi secondi messi in fila risulta 
chiaro chi vive e chi muore tra il leone Ln, e il me stesso 
Qn, e nei secondi seguenti stanno certamente svolgendosi 
[...]. Ogni secondo e defmitivo, chiuso, senza interferenze 
con gli altri. 



I. Calvino, Ti 



con zero 



e will now introduce a completely different approach to the 
tomography problem. Let us consider the image recon- 
struction as a problem of interpolation. In fact, given a 
FfiST ^y basis of functions {&«(a;)}j=i...„ that interpolate the func- 
tion /, i.e. such that 




f( X ) = ^2°i b i 



X) 



Vxel 



i=0 



where X is properly chosen, then for linearity of the Radon transform shown 
in Lemma (3.1.2) 

n 

Kf(y) = $>7^(y) Vy G Y = {(t 3 ,9j)}. 

i=0 
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Now consider the problem from a practical point of view: incoming data are 
stored in a matrix / x p where / is the number of angles of acquisition and p 
the number of pixels recorded for each angle ($i, . . . ,$i). We want to obtain 
as output an image oi N x N pixels. Let us consider now the domain as this 
grid of pixels: let Pi = [xi,Xi+i] x [j/i, 2/i+i] the i-th pixel, where (xi,yi) is its 
top-left vertex. Then we introduce the following basis 

h(x,y) =X Pi ( x iV)- 
We know the Radon transform of each one of these items 

nbi(t,6) = meas(£ tfi D Pi) 

where £ t> g = {t8 + s9 L \ s G M}; hence this integral can explicitly calculated 
for every angle chosen. 
The obtained system is 

Ax = c (4.1) 

where A(i, j) = lZbi(tj,0j) is a sparse matrix N 2 x Ip, x is the unknown vector 
such that Xi = f(xi) and c is the vector of the projection data Cj = lZf(tj, 9j). 
The methods using this approach are known as Algebraic Reconstruction 
Techniques or ART. 

4.1 Some iterative methods 

Nowadays the obtained system is solved numerically in a large amount of 
time: in fact if the matrix has standard dimensions (TV = 128, I = 60, 
p = 128) and to store it we need more than 1Gb of memory. The problem 
can be bypassed considering the sparse structure of the matrix A (at least 
90 % of its elements is equal to 0). For solving the system (4.1) we also will 
use iterative algorithms, whose general structure is the following: 

• the initial vector x^ is a blank (the zero vector or unit vector) or a 
inaccurate reconstruction. 

• the image at step k, x^ k \ is projected and compared with the data (i.e., 
for instance, e^ = c — Ax^). 

• the image is modified considering the error found in the previous step 
^(fc+i) _ g(x( k \eW) with g a certain function). 

Here are some of the most common algorithms: 

32 




Figure 4.1: The main idea behind Kaczmarz algorithm is the projection on the 
columns of the matrix A; in the figure we used three columns as example. 

4.1.1 The Kaczmarz algorithm 

Let us choose iteratively a row k of the matrix A (usually % = k mod Ip or 
we choose randomly). Then we project the image on the vector A^. 



x 



(fc+i) 



X 



(k) 



\T 



-(*) 



+ 






In practice using this algorithm we reach a satisfying reconstruction only 
after having projected the image on all the rows of the matrix. In standard 
dimensions, the solution is accurate only after O(10 4 ) iteration. That is why, 
in order to increase the speed convergence, sometimes we adopt a relaxation 
parameter A, chosen properly or randomly (or with probability proportional 
to II A- II 2 ) and the iterative scheme becomes 



x 



(fc+i) 



x 



(k) 



+ A* 



Ai 



x 



(k) 



AT ■ A,. 



-Ai 



in this case the method is called Kaczmarz algorithm. It is shown in [18], 
[13] Ch 5.3 pl30-132 and [10] Ch 2 p. 2 that the convergence is guaranteed 
for Afc G (0.2). The ART algorithm is the first method historically used, but 
it was abandoned for the FBP because of its slowness of convergence. 
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4.1.2 Maximum Likelihood Expectation Maximization 

The main idea of this algorithm 1 , based on probabilistic arguments, is to find 
the x that maximizes the likelihood of the probability L(x) = P(c\x). This 
means that we are looking for the solution x which maximizes the conditioned 
probability of having our data c with an image x. 
Assuming that the noise is Poissonian on c, then the likelihood is 

lP «) C \- C? 



L x) = P ex) 



n 



Oil 



where c* = Ax is the exact sinogram, i.e. the projection of the exact solution 
x. Equivalently, it maximizes 

i P 
l(x) = log(L(x)) = J2 ~( Ax )i + °i log(Ar)« + K 



i=l 



with K a constant and constraint no n- negativity on X{. Kuhn- Tucker condi- 
tions are then: 

dl(x) 



3 dx,; 



dl(x) 



dxi 



if xj > 



< if xj = 



so that we get from the first equation 

/ 



dl(x) 

= *1Z7 



Xj 



\ 



l P 



lp 



-J2 ai, j + Yl 



QijCi 



i'=l 



i=l 



N 2 

E 



tX/r} n ' tXj -1 ' 



or equivalently 









ip 



J2 ar 



i P 
E 

2=1 



^i,j^i 



N 2 



■' -J 



from which comes out the iterative scheme 

M i P 



j/Xji 



x (k+i) = x j y^ 



(X{jC{ 



ip 



N 2 



i'=i j'=i 

1 see [18] for a more detailed explanation. 



W 
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Iterative reconstruction after k^l iterations 



Iterative reconstruction alter k=3 iterations 



Iterative reconstruction after k=5 iterations 






iterative reconstruction after k=1G iterations 



Iterative reconstruction ^ttci k=20 iterations 



Iterative reconstruction after k=50 !tvr3: : 0!:s 






Figure 4.2: Several iterations in MLEM reconstruction of a phantom. 

We now rewrite this scheme in a vectorial and multi-step form 

1. c f := AxW 

2. c q := c./c-f (punctual division) 

3. x b = A T c q 

4- S j = l^i a i,j 

5. x( k+1 ^ = x( k \ * x b ./s (product and division are made elementwise) 

Faster results are obtained by using the OSEM algorithm (Ordered Subset 
Expectation Maximization); we simply apply EM reconstruction of ordered 
subset of all the projections {9i}i=i...i — © = (J Qj. 

The OSEM algorithm is currently the most used in nuclear medicine 
for its speed of convergence and accuracy. 
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4.1.3 Conjugate gradient 

The conjugate gradient method is widely used in all fields for solving linear 
systems Ax = b with A square matrix of order n, symmetric and positive 
definite. In our framework the method is known as LSCG (Least Squares 
Conjugate Gradient) because we apply the conjugate gradient method to 
the normal equation of the system (4.1). 

In fact the matrix A T A is square, symmetric and positive definite. The 
method proposed below solves iteratively the normal equations 

A T Ax = A T c 

without having to compute explicitly A T A. The algorithm starts with an 
initialization phase: 



• given x^ ' initial value 
. s (o) = c _ Ax (o) 

m r (o) = p (o) = at s (o) 

• g(°) = V°> 

then begins to calculate, at each step: 

| r M|2 

• r (k+l) = J\T s {k) 



P 



| r (fc+l)|2 



| r (fc)|2 

p(k+l) _ r (fc+l) _|_ gp(k) 



q 



(fc+i) 



Ap (k+i) 



This algorithm, despite the appearance, is in most cases as fast as an 
iteration of EM. 
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4.2 Accelerated reconstruction 

Although iterative algorithms are extremely fast, it is possible to further 
accelerate the reconstruction, by using a technique of zero-padding of matrix 
A, i.e. 

_ / a i ,3 if a i,j > 7 

a{l ' j) ~ \ else 

The algorithm is trivial, 

A = A. * (A > 7) 

where .* is again the pointwise multiplication. This algorithm sets to zero 
all values of A under 7, increasing the sparseness of the matrix and conse- 
quently the velocity of resolution of the system (4.1); it is evident that this 
reconstruction is faster but less accurate. 

To choose the parameter 7 we can proceed in different ways: we can put 
7 = e, the machine precision, if we want to delete only the items nearest to 
zero and not decrease the accuracy of the solution. More typically we choose 
a fraction of the maximum value of A such as 7 = 0.05 max{ojj}. 
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Chapter 5 
Kernel methods 




n alternative choice of the basis in the interpolation can lead to 
different results in terms of time of computation and accuracy 
)t?) of the solution. Let us recall the interpolation method. Given 
a basis of functions {bi(x)}i=i... n that interpolate the function 
/, i.e. such that 

n 
i=0 

then for the linearity of the Radon transform 

n 
i=0 

We will see in this chapter the interpolation problem solved through the use 
of kernel functions as basis. 



5.1 Direct kernel interpolation 

To introduce the new basis we need some definitions. 

Definition 8. A function f : Q C ]R S — y M. is called essentially compact 

supported if the closure of the set S — {x s.t. \B(x)\ > e} is compact 1 in 

R s . 

Remark: if we choose e = the definition above is the classical definition 
of compact supported function. 



with e chosen properly, e.g. the machine precision. 
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Figure 5.1: Original phantom and kernel reconstruction with Gaussian kernel and 
shape parameter e = 1 after 50 iterations of LSCG. 



Definition 9. A function f : Q C ]R S — y R is called radial (in this work we 
will also call them Radial Basis Functions, briefly or RBF) if 



<t>{x) = tp{\\x\\) 



where tp : [0, +oo] — ¥ R. 
A function K : Q x Q C 



— >■ 



zs called Kernel if 



K(x,y) = <j)(x-y) 



with (j) a radial function onnCR s . 

Instead of the natural pixel basis we use another basis of functions bi(x) = 
B(x — x~i), where B is essentially compact supported and radial. For Lemma 
(3.1.4) the Radon transform of a radial function does not depend from the 
angular variable 6, so to build the matrix we only have to compute lZB(tj) 
where tj = d(£f,e ,%i), with d the euclidean distance. Since the line t t p = 
(x(s)) = (tO + s^- 1 ) has the equivalent, nonparametric form t = xcosO + 
ysin.6, so the distance from the point x = (x,y) to the line £ tj g will be 



U 



d(i th e p Xi 



|xcos^ + ysin.9 — t\ 



In Table 5.1 we show the transform of some RBFs, calculated in order to 
use them in interpolation. Since the matrix of the interpolation system will 
be very large, we have chosen (essentially) compact support functions, so 
that the system matrix will be sparse and the solution of the system will 
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be computed using one of the iterative algorithms described in the previous 
chapter. 



Function name 


f 


Kf 


Ball 


Xi?(i,o)( r ) 


y/G- t2 h 


Gaussian 


2 2 

e- £r 


Vl e -e 2 t 2 

E 


Wendland y? 2 ,o 


(l-er)l 


{2eV^ 2 + §eV 3 " .log(2^)t 2 }X h£ , £l 


Wu ipi,i 


(1 -er)\(er + 2) 


|.3iog( 2 ^)tnx h£ , £l " 



Table 5.1: Radon trans form of the four chosen functions. 

The shorthand y/~- = \/(-s — i 2 )+ so that t G (0, -) and r = \\x\\, as usual. 



5.2 Definite positivity 

In the theory of interpolation with kernel functions, which can be found in 
[4], a typical request is made on the chosen kernel, made to ensure that the 
interpolation matrix will be nonsingular and positive definite. This condition 
is just called positive definiteness. 

Definition 10. A symmetric kernel function K is called positive definite 

if\/X — xi, . . . , Xn C Q data points set and c G M. N 

N 



^ c i c j K(x i ,x j ) > 



and it's called strictly positive definite if the inequality above is verified 
and the equality stands only for c — 0. 

It is immediate to observe that the matrix of interpolation of a strictly 
definite positive kernel is definite positive (semi-definite if the kernel is only 
definite positive). 

In our setting the system matrix is not K(xi, Xj), but TZ(K(x, Xj))( yi y So we 
have to verify that the Radon transform of the chosen function is positive 
definite. Our next step is a theorem that proves that the Radon transform 
preserves the definite positivity of a radial function. The next theorem will 
be useful for this purpose. 
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Theorem 5.2.1. Suppose G L 1 (IR d ) and a continuous function. Then <f) is 
positive definite on IR d if and only if is bounded and its d- dimensional Fourier 
transform J-d(4>) is nonnegative and not identically zero. 

Proof, see [15] Ch 3 for further details on this theorem. □ 

Using this result we prove the following Theorem, which is fundamental 
in our applications. 

Theorem 5.2.2. If (f>{x — y) = K(x,y) is a radial function, <fi G L 1 (IR (i ) ; 
continuous, bounded and positive definite on IR 2 ; then its Radon transform 
TZf(t) is bounded and positive definite on M. 1 , provided IZf G L x { 



Proof. A direct consequence of Theorem (5.2.1) is that ^(0) (P) is nonnega- 
tive and not vanishing, and this property remains if we consider the function 
in polar coordinates ^{^{T cos 6, T, sin 8). By means of the Central Slice 
Theorem (3.3.2), also the function 

Ji (ft/) (T, 0) = .F 2 (0)(Tcos0,Tsin0) 

is nonnegative and not vanishing. Now we remember that the Radon trans- 
form of a radial function depends only on the linear variable t, then 

T\ (ft/) (T) is nonnegative and not vanishing 
and, again using Theorem (5.2.1) we have the thesis. □ 

All functions we are using are positive definite on IR 2 , and so will be their 
respective transforms. The only exception is the ball function, which is only 
conditionally positive definite 2 of order 1; we have included that function 
only to compare the results with other kernels, but we have no warranty 
about the nonsingularity of the system matrix. 

5.3 Error estimate 

We want to estimate the reconstruction error in case of direct kernel inter- 
polation; if we call such error 



e(x) = f(x) - y^ CjK, 



.X) 



2 a weaker condition, see [4] for the definition. 
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then for linearity of the Radon transform 

e(x) = IT 1 (lie) (x) = IT 1 (llf(y) -J^c^K^y)) (x) 
and so 

n 

IMIl 00 ^ 2 ) < ||^ ||(oo,oo) \y^f ~ / J c t'^-^"ilU° c (Rx5 1 ) 

i=0 

where, as in the matrix case, 

||^ ||(oo,oo) = SUPg^L * 



\\9\\oo 

Now, if we want to proceed with the error estimate we need to show some 
elements of the approximation theory in interpolation with kernel functions. 
For a more detailed panoramic about RBF interpolation, see [4] in particular 
Ch 14. 

Definition 11. Given a set X of N pairwise distinct points of Q C 1R S we 
call fill- distance the radius of the largest ball in Q that does not contain any 
point of X, i.e. 

h = hx,n = sup min \\x — Xj\\2 

Definition 12. Let H be an Hilbert space H — {f : Q C M s — > M.} with 
scalar product (•, •). A Kernel K : f2xf2 — > R is called reproducing kernel 
for H if: 

1. K(-,x) EH ViGH 

2. f(x) = (f,K(.,x)) \/feH,xen 

Definition 13. Given K radial basis function on IR S strictly positive definite 
its native space is the Hilbert space Nk{&) = Hk(VL), where Hk{Q) = 
span{K(-,y) | y G Q} is a pre-Hilbert space with K its reproducing kernel. 
Its scalar product is the bilinear form 



N N \ N 

/(•) = ^2ciK(-,Xi),g(-) = ^2jiK(-,Xj) ) =Y1 c i c j K ( x h 

i=\ jr=l / R i,j=l 

given a set X of N pairwise distinct points ofQ. 
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Xj) 



Theorem 5.3.1. Let Q C IR S be an open, bounded set which satisfies the 
interior cone condition 3 .Suppose K e C 2k (Q xQ) is a symmetric and strictly 
positive definite kernel. Denote the interpolant to f E jVr-(JI) on the set X 
of N pairwise distinct points of Q as Vf. Then there exist positive constants 
ho and C such that, if hx,a < ho, then 

\f(x)-V f (x)\ < Ch x ,nVQM\\f\W Km 

where C k (x) = max^^ raax W:ZennB{S:C2hxn) \D^K(w,z)\. 

Corollary 5.3.2. With the same hypothesis as above, there exist positive 
constants ho and C such that, if hx,n < ho, then 

Wf-VfWoo^Chx^WflW^n) 

Now we can formulate the theorem that will give an estimate of the error 
on the image reconstruction problem. 

Theorem 5.3.3. Let f G C^°(B (0,1)) a b-band-limited function, and let 
g = IZf be reliably sampled. Let K be the interpolating Kernel function such 
that 1ZK and its domain satisfy all the conditions of the Theorem (5.3.1). 
// we call the interpolation error 

n 

e(x) = f(x) - ^ c i K i{%) 

i=0 

then there exist positive constants ho and C such that, if hx,n < ho, then 
1 1 e 1 1 z,oo < 2\S | ||wjrf ea z||x,i Chx,n\\ T^f\\M nK (Ci) 

where \\w idea i\\ L i < &yf§. 
Proof. As we have seen before 

n 

||e||L°°(K 2 ) < 11"^- ||(oo,oo) WR-f — y jCjIZKjWlcg^s 1 ) 

where the first part of the second term is 

ii-o-in l|7^~ lo lloo 

\\7t (00,00) — SUPg^Loc - - 

\\ n\\ 

\\y iioo 

3 or more simply dfl has regularity at least C 1 . 
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and the second one is simply an error of the interpolation using a positive 

definite kernel. 

So all we only have to estimate the term ||7£ -1 || using known results for 

analytical reconstruction. 

From Theorem (3.4.3) we get 

||/||l°°(r 2 ) < 2\S I ||wft|| i i( R )||g||Loo( Rx>g i) + |e 3 |. 

In our formula the transform of the kernel function is made analytically, so 
the term e% is null, and instead of ||iUb||i,i we have to use the ideal filter 
w = Wideai such that lZ*w = S, the Dirac function. Let's think about the 
known inversion formulas for the Radon transform. 



f = n* 



r-'C-ffW)) 



K* [w * Kf] 



Let us call A = w * TZf; then the Fourier transform gives 

F(A) = T(w)T(Kf) 
from one side. On the other side 

.F(A) = ^{JZf) 
So, in conclusion, 

The second term would be divergent, but remembering that we are consider- 
ing a 6-band-limited function, so that the filter must be limited by the radius 
of acquisition r and by b in the frequency domain 4 

IM|i.i(R) = IMU 1 ([-»-,r]) < v / 2r||w|| L 2 ([ _ T , ir]) 

using the Parseval's identity on the equivalent form \\gW2 = -k=\\J r ~ 1 g\\2 and 
knowing that the radius of acquisition 5 r = iV^ with iV as usual the number 
of pixel per dimension we get 



2irN / 1 \u\ N 2b 3 N 

W\»w < V ~w ' V 2^ " Y h2 ^ b ' b]) = \lm-\lY = b \l 18 

the inequality is obtained using a known embedding theorem of L p spaces on limited 
domain, if < p < q < +00 and the domain of / E LP is a limited set S \\f\\ p < 
meas(S)p~* \\f\\ q . 

5 remember that the size of a pixel is FW ^ M = ||. 
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This means that 

'N 




II^IUoc) =21^1 \\w\\ L i(s ) <2\S 1 \b\ 

In conclusion 

n 

||e||z,°°(IR 2 ) < ll"^- II (00,00) WR-f — / ,CjR'Ki\\L°°(l&.xS 1 ) < 

i=0 

< 2\S I H^idea/IU 1 Chx,n\\Tlf\\Af nK (n) 
which is the promised result. □ 
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Chapter 6 

The emission tomography case 



Denise and Steffie stayed home that week as men in Mylex 
suits and respirator masks made systematic sweeps of the 
building with infrared detecting and measuring equip- 
ment. Because Mylex is itself a suspect material, the 
results tended to be ambiguous and a second round of 
more rigorous detection had to be scheduled. 

D. DeLillo, White Noise 

s we have said in the previous chapters, for reconstructing 
an image in the case of emission tomography we will have to 
deal with attenuation factors. Two different approaches can 
be used. The first one is to reconstruct the image with the 
method known in CT, using an attenuation correction method (more used 
in the PET), with the help of some algorithm for the segmentation and the 
detection of the contour of the region of interest (ROI). 
The second consists of inverting the attenuated Radon transform once we 
have estimated the attenuation coefficients a(x, y) from a SPECT/CT hybrid 
machine (or possibly SPECT/MRI 1 , with greater difficulty). In this chapter 
we will discuss about the attenuation correction methods and the analytical 
and iterative algorithms for the inversion of the attenuated Radon transform. 

6.1 Geometric methods for attenuation correc- 
tion 

The resolution of the reconstruction problem in the case of emission tomog- 
raphy type is complicated because we do not really know the attenuation 




1 MRI, Magnetic resonance imaging, another technique used in radiology. 

47 



coefficients of the different tissues within the body. To overcome this prob- 
lem we use two algorithms of attenuation correction which are based on a 
geometrical estimate of the coefficients. 

6.1.1 Soreson algorithm 

Soreson proposed a pre-correction of the data obtained from the projections. 
The data is in fact multiplied by the correction factors in order to obtain for 
each angle the number of counts that would be obtained without attenuation. 
It is hypothesized the uniform distribution of the tracer and the uniform 
density of the tissue along each ray; this means that the algorithm will give 
out good results only when the area to reconstruct is sufficiently uniform, 
for example the brain. Let P^q and Pk,e+n two opposite projections, T the 
diameter body and t the diameter of the active region along the straight line 
projection. 



Then, if we call N = yPk,e • Pk,e+w, the geometric mean of the data the 
correct value for both angles is 



sinh(/7 T 



-.W Lll 



where 7 is the linear attenuation coefficient (e.g. the brain, if the tracer used 
and 99m T we suppose 7 = 0, 035cm~ 1 ) and / = ^ is a factor which in practice 
varies from 0.5 to 0.75. We can fix that value, for instance to 0.6 or better 
to an average estimated value of ^. 

6.1.2 Chang algorithm 

This algorithm makes a post-correction on the reconstructed image. The 

result of the classic reconstruction is multiplied pixel by pixel by the weights 

obtained from an estimate of the mean attenuation that the beam with angle 

9 had been subject to at the transition from pixel (x,y). 

Then, if P is the image obtained e.g. by filtered back projection and Pq the 

corrected image, then Po(x,y) = c(x,y)P(x,y), where the correction terms 

are: 

C(*>y)=(yE' 

where lg i is the distance between the pixel (x, y) and the edge of the body 
at the angle 9i and /i is the linear attenuation coefficient, supposed to be 
constant. 
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1 ■ - l 



N 



T 



} 



•> 




P(t,n + 6) 
Figure 6.1: Scheme of the data used by the algorithm of Soreson. 



Therefore also the Chang algorithm gives good results if the tested area is 
uniform 2 . 



6.2 Inversion formula of the attenuated Radon 
transform 



We recall the AtRT formula 

K*f{t,6) 



-Va(x,0+l) 



m 



(6.1) 



(t,9) 



now, given a(x) we want to recover the value of / from its transform g = lZ a f . 
The first analytical solution to this problem was given in [14] by G. Novikov 
in 2000. We will see the version of this formula published in [12] by F. 
Natterer in the same year. 
First of all we need the definition the Hilbert transform. 

Definition 14. Let g a suitable function, then its Hilbert transform is the 

function 



Ug{s) = l[M df 



7T Jr S - t 

where the integral is meant as a Cauchy principal value. 
2 for a more detailed explanation, see the original article [1]. 



(6.2) 
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Note that this transform can be thought as the convolution g * h where 
h(s) = — , but the convolution integral does not converge. 
Now let us define the function 

h=-{I + iU)na (6.3) 

where I is the identity operator, i the imaginary unit, and both the operator 
are meant to be applied on the linear variable of IZa. 
Now we can enunciate following Theorem 

Theorem 6.2.1. (Novikov-Natterer formula) Let g = lZ a f for f trans- 
formable function, and h the function (6.3). Assume a(x,y) known, then f 
is uniquely determined by the following formula 

f(x) = ^t div J ^ Qe v ^^\e- h Ue h g)^. m dQ (6.4) 

where S x = [0, 2tt]. 

Proof, see [12] for the complete proof. □ 

Now let a = everywhere, then g = IZf, h = and T>a = so that 

which is an equivalent form 3 of (3.3). 



6.3 Discretization of the AtRT inversion for- 
mula 

The discretization 4 of the (6.4) is similar to the discretization of the FPB 
formula (3.4). In fact all we have to do is to take the "attenuated back- 
projection" of the convolution product between g and the inverse Fourier 
transform of a low-pass filter. 
First of all we have to compute the function g a 

g a = y\z{e- h Ue h g) 

Putting h = hi + ih,2, h± = \lZa and hz = \T-L1Za then 

g a = e~ hl (cos hiHe hl cos h 2 + sin h 2 l-Le hl sin h 2 )g 



3 the proof is in [5] Ch 6 p. 51. 

4 see [12] Ch 3, p. 6 for more details about the discretization. 
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where, like in the non-attenuated case, % is approximated by a convolution. 
Let o"6 a low pass filter 5 with b its cut-off frequency and £5 = J 7-1 ^, then 

H b f = S 6 * / 
once we have calculated g a we only have to compute 

f(x) = ^div J ^ 6e Va ^ e+ ^g a (x ■ 6,6) d6 = n* a (g a ), 

where the integral is computed in the sense of a discrete mean and the di- 
vergence operator is approximated by finite differences as explained in [12]. 



Attenuation phantom 



Activity phantom 





Reconstruction of the attenuation map 



Reconstruction of the activity 





Figure 6.2: Analytical reconstruction of a SPECT/CT phantom data. 



6.4 Iterative methods 

We propose a new method of reconstruction for the attenuated case. The 
main idea is the use of iterative methods as in the transmission tomography. 

5 again, for a list of the low-pass filter used see the chapter about regularization. 
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i=0 



We consider the outgoing Gamma-rays at angle 6 coming from a pixel with 
unit concentration of tracer. 

As we said in previous chapters, the AtRT is not equivalent to the RT of a 
weighted function. Nevertheless from the matrix A used in the transmission 
case, we can obtain a matrix B such that Be = d where Bij = lZ a bi(tj, 9j), 
Ci = f(xi) the unknown term and dj = TZ a f(tj,9j) the data. 

Let us consider, as before, the basis bi(x) = X p and assume 

n 

/(£) = 53 Cibi 

Then for linearity 

n 

K a f{y) = J2 CiHabiiy) VyeY = {(£„ 0,)} 

i=0 

Now to compute the elements of B, we have to remind that the data of 
attenuation come from a CT reconstruction, i.e., as suggested in [7], are on 
the form 

N 2 

fc=l 
For linearity its fan beam transform will be 

N 2 ^+oo ^ 2 

Va(x, 6) = J29k X Pk (s + ^) dt = J29k meas(P k n £+ e ) 

*=i Jo fc=i 

where £f e — {x + W \ s E 1R + } and meas is the Lebesgue measure. Now let 
us consider the AtRT of the basis hi 

K a bi = exp I - y^ g k meas(P k n ^Je+f ) ) Xp,- 

^^.» \ fe=l / 

This last formula is hard to compute, but it offers a useful interpretation. In 
fact if the tracer is concentrated only on the i-th pixel Pi, the outgoing rays 
will be subjected to an attenuation equal to the attenuated sum of the path 
length; if I in — 1, according to Beer-Lambert's law 

/ N 2 \ 

hut = Iin exp I ~^2 9k meas(P k n £± ) J . 
52 



Now, if we consider the matrix A used in CT tomography, we can compute 
the outgoing rays from the pixel P t in (tj, 6j) as 

Bij = A itj exp ( - ^ 9k meas(P k r\V. tjfij 

Aij exp | ^2 9k A klM 

(fci,fc 2 )e^(ij) 



where K, 



(id) 



{(fci, fe)} C {1, . . . , N } is the set of the indexes of the pixels 



covered by the line 



and the index k = Ip — {{k\ — l)p + fo) + 1- 



This method requires a computation time a little longer than that of the 



Attenuation phantom 



Activity phantom 





Reconstruction of the attenuation map 



Reconstruction of the activity 





Figure 6.3: Iterative reconstruction of a SPECT/CT phantom data with A = 0.1. 

analytical approach because the system matrix B is calculated a posteriori. 
Moreover this matrix has often a large condition number, so the analytical 
solution seems to be preferable this time. 

Nevertheless, a regularization of the linear system can bring to a more accu- 
rate solution than the analytically reconstructed image. First of all we can 
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introduce a relaxation parameter A G [0, 1] in order to consider relatively the 
effect of the attenuation 



B iJ = Ai J ex P \~ X 5Z 9k Ak ^ k ^ 

Note that B^ = A and B^ = B. With this little foresight the linear system 

can lead to a more accurate solution. 

Regularization methods will be introduced in the next chapter. 
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Chapter 7 

Regularization methods and 
post-processing 



)^"?W™^(1 he matrix of the linear systems to solve both in transmission 
c^\ B/^ an< ^ em i ss i° n cases is large, sparse, and often ill-conditioned. It 
^^BfxS' can even happen that the matrix had pair of identical columns 
^tIMI^^ and some null columns. Moreover, the scattering can bring to 
a noisy data. So it is reasonable to add a regularization in the reconstruction 
which will bring an image with less artifacts 1 . We also introduce a contour- 
detection algorithm. 

7.1 Low-pass filters 

Let L be the cut-off frequency e.g. the resolution frequency b found in equa- 
tion (3.2). Now we will choose the filter wl, and convolve its Fourier trans- 
form with the reconstructed image P. The result will be a filtered version of 
P. 

P = T- 1 [w L T(P)] = P*T~ 1 w L 

The following functions are commonly used in filtering medical images: 
• "Ram-Lack" filter 

^(IHI) = X[-l,l](IHI) 



1 a regularization method in general increases the smoothness of the solution, i.e. in 
our case increases the partial volume effect and the number of artifacts can increase. So 
wee will always use regularization methods in a further step, instead of embedding it into 
the reconstruction process. With this expedient we will have two images to compare: an 
unfiltered, raw image with better resolution and a filtered, smooth image, but with a lower 
resolution. 
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which is simply the truncation of frequencies larger or smaller than L. 
"Shepp-Logan" filter 




M\\M\)= MkSK Xi-^i(HHII) 



2L 



essentially a truncated sine function, 
"low-pass cosine" filter 



^(IIHII)= COS (^l) X[-L,L](IIHII) 



a normalized, truncated cosine function. 

In our notation we wrote |||^||| not as a norm but as |||^||| = max{i/i,^}. 
The functions used as filters in FBP formula and its analogue in the emission 
case are exactly the version of those filters in R 1 (then we use v instead of 
lll^lll) multiplied by the absolute value of v. 

u L {v) = \v\w L (v) 

If we want to find the discrete value of the cutoff frequency L we have to 
know the following proportion L : f s — cut : N where cut is the discrete 
cutoff, N the numbers of samples and f s the sampling frequency, which for 
a standard dimension problem is f s = l/3.5mm _1 = 0.286mm _1 because 3.5 
mm is the typical value for the physical length of a pixel 2 . 
The results of filtering in our experiments are not encouraging. In fact filter- 
ing does not bring to a better image, either in terms of oo-norm error or in 
terms of visual quality. 

7.2 Philips-Tichonov regularization 

Philips-Tichonov is a classic regularization method used for ill-conditioned 
system, its main purpose is to add regularity to the solution; the typical 
iterative algorithm for solving linear systems minimizes the quantity 

min 1 1 Ax — 6IU 

x£R n 



2 more precisely the length of a pixel is FWHM/3, and the typical value for FWHM is 
about 10mm around 150mm of distance. 
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now we add a penalty term 

min \\Ax — 6II2 + 'Y 1 1 -Z^a; 1 1 2 

xeR n 

where 7 is a positive, tunable constant, and F = A m the regularization ma- 
trix that makes the solution similar to a function of regularity C m . Typically 
Aj.j = 1 and A,; ., = — 1 if the i-th pixel is next to the j'-th, zero otherwise 3 . 
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The solution of the new system will be 

(A 1 A + 1 F t F)x = A% 
and we will find this solution iteratively. 



Iterative reconstruction af.ei lo-100 iterations 



sn and reconstruction after k=1 00 iterations 





Figure 7.1: Reconstructed image (left) and its regularization (right). As we can see 
from the oo-norm errors a regularization does not always bring a great improvement 
to the quality of the solution. 



7.3 A segmentation algorithm 

A phantom is a function to test the efficiency of the algorithm in silico 4 , 
i.e. without data given by a physical machine. Is is also called phantom a 
physical device that simulates a certain organ and used to make experiments 
on the machine. For clarity we will call the second one dummy. 



3 we assume A = I, the identity matrix. 

4 this term is uses commonly in biomedical engineering; it means that experiments are 
not made on a living been (in vivo) or on a few component or chemical reproductions (in 
vitro), but only on a mathematical model. 
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Usually a phantom is a step function, or rather / = ^V \^ B _ with Bi C R 2 
some compact sets. This type of functions is congenial to the reconstruction 
because it splits the region of interest (or ROI) in different areas where the 
attenuation coefficient is constant. We can assign constant that if we recog- 
nize the type of tissue or calculate in the case of SPECT/CT. Therefore we 
need a segmentation algorithm that levels the image obtained by the CT im- 
age reconstruction, making it like a phantom, which we call phantomization 
algorithm. The proposed algorithm is the following: 



exact phantom 



EM reconstruction 



"phantomized" reconstruction 






2-norm error: 36.6461 



2-norm error: 0.10308 



Figure 7.2: Reconstruction of a David star phantom and its phantomization; as we 
can see from the errors the phantomized solution is close to a phantom. 



1. P is the NxN matrix of reconstructed image, k the number of different 
levels that you want to look in [0.1] 

2. put its elements in a column vector, P = P(:) long N 2 

3. for i = l,...,7V 2 -2 

(a) we calculate S the sum and D the difference of the elements in 

X = P(i:i + 2) 

(b) if S or D has norm less than | then P{i : i + 2) = mean(X) were 
mean is the the average 

(c) otherwise P remains the same 

(d) end of the for loop 

4. repeat 2. and 3. for each rotation of | of the matrix P 



58 



This algorithm, if the starting reconstruction is sufficiently accurate, returns 
an image very similar to a phantom. We can alternatively assign to P dif- 
ferent values as P(i : i + 2) = 1/k * roundik * mean(X)) where round is a 
rounding function. In this way the zones are more uniform, but the value of 
the attenuation coefficient tends to be less accurate unless we choose k ad 
hoc for the values of attenuation coefficients in the examined area; e.g. if 
we know that the coefficients are all integral multiples of ^. For instance 
k = 10 is a good value. 



7.4 A simple algorithm for the contour detec- 
tion of the ROI 



sinogram 




iterative reconstruction at height=92 on 1 23. 



contour of the ROI 





Figure 7.3: Edge detection of a brain dummy's data; the algorithm detects effi- 
ciently the ROFs contour. 

There are many algorithms written for the detection of the contours of 
an image, the one we are proposing is a very simple application of the convex 
hull algorithm of a point set. 

Let x be the solution e.g. from an iterative reconstruction. Now we identify 
the area in which is the traced is more concentrated with a simple zero- 
padding the solution x = x. * (x > 5) where S = 7 max(x) and 7 G (0, 1) is 
chosen properly (e.g. we have chosen 7 = 0.22 in our experiments). Now we 
apply a convex hull algorithm 'convhull' and draw the segments between the 
external points. 
The following MATLAB® function builds a matrix B with all zeros except 
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for the contour points, where the matrix has value 1. 

function B=edge_detection(I ,rapp) 

°/ohome-made edge detection for SPECT reconstructed images 

B=zeros(128,128) ; 
m=max(I)*rapp; 
II=zeropadding(I,m) ; 
II=reshape(II,128, [] ) ; 
[x,y]=find(II>=m) ; 
if length(x)>2 
k=convhull(x,y) ; 

lung=30; %lung is chosen properly to represent every pixel of the line 
t=linspace (0 , 1 , lung) ; 
for i=2: length (k) 

a=round(t.*x(k(i-l))+(l-t) .*x(k(i))) ; 

b=round(t.*y(k(i-l))+(l-t) .*y(k(i))) ; 

for j=l:lung 

B(a(j),b(j))=l; 

end; 
end; 
end; 
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Chapter 8 

Experimental results from in 
silico tests 



When I was a child, I used to talk as a child, think as 
a child, reason as a child; when I became a man, I put 
aside childish things. 

At present we see indistinctly, as in a mirror, but then 
face to face. At present I know partially; then I shall 
know fully, as I am fully known. 

ICor, 13: 11-12 




e now present some experimental results from several tests 
made on phantoms (in silico). The aim of these experi- 
ments is to test the algorithms presented in this work and 
to compare the quality of the respective reconstructed im- 
ages. All the tests has been made on the standard dimensions, i.e. we have 
reconstructed a 128x128 Shepp-Logan phantom by taking 3°-step angular 
projections of 128 linear data each, except where is differently specified, and 
the error will be evaluated in oo-norm, coherently with the error estimates 
made in previous chapters. 

The computer used for all these tests have a 3.5 Gb of RAM equipped with 
two dual-processors of 2.70 GHz each on a 32-bit architecture. 



8.1 Filtered Back Projection results 

We now proceed by evaluating the performs of the FBP algorithm in terms 
of precision and time of computation. In these experiments we are using 
the standard radon and iradon functions included in the image processing 
toolbox package of MATLAB®, which standards are different: it reconstructs 
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128x128 images by a 185x60 matrix data. Now, to test the error bound 
provided by Theorem (3.4.3), which is 

11/ — /||l°°(k 2 ) < 215 | ||w&||i,i(R)||<7 — (?||l°°(RxS' 1 ) + |e 3 j , 
we have to compute the L 1 -norm of the filters in use. 

• "Ram-Lack" filter 

\\w L (v)\\ L i =2 [ v dv = L 2 = c (1) L 2 
Jo 

• "Shepp-Logan" filter 



L / 7r||i/||\ , 8 



o 



w L {v)\\ L i = 2 / v sine — — dv = —L = c {2) L 



2L 



2 „ r2 



TV 



• "low-pass cosine" filter 

\\wM\W = 2 J L ucos [^j dv = fc^L 2 = c (3 )L 2 

All the constants cu\ are 0(1), so the different choice of filters does not 

change the order of the error bound. We choose the constant C3 = 0.3974 as 

the FBP error with no noise on the signal. 

In Figure 8.1 is shown the error of the FBP algorithm using the "Ram-Lack" 

filter for several noise errors. 

The FBP algorithm is quite fast, its computational cost is 0(iV 2 log iV). 

The mean time of computation for this experiment is t = 0.088s. As we can 

see from Figure 8.2 the FBP can reconstruct quite large images in less than 

a second. 
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Figure 8.1: A test for the error bound using the standard "Ram-Lack" filter; note 
that the unfiltered back projection error is bigger than the bound for low noise 
levels, and reaches the same order of the bound for higher noise levels. 



computing 
time of 
FBP 
(s) 




Figure 8.2: Time test for the FBP algorithm. 
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8.2 Iterative Methods results 

Next experiments concern the iterative methods described in Chapter 4. 
Since the ART algorithm is very slow 1 we proceed by showing the results 
of MLEM and LSCG methods. The algorithms in these experiments are 
exactly the ones described in the corresponding, sections, while the algo- 
rithm used for computing the system matrix A is contained in the "Air tools 
package", [8] available in Hansen's web page 2 . 

In this section, we will show some tests on several iteration numbers and 
noise levels. 

In the case the noise on the data are null, both algorithms show fast and 
good performances. The results of our test are shown in Figure 8.3. In both 
cases the time of computation grows linearly with the number of iterations 
and the error is almost everywhere decreasing. This confirms that the com- 
putational cost for k iterations is 0(kN 2 ). 

We note that MLEM stops after the 75th iteration. In practice, MLEM al- 
gorithm is more precise than the LSCG, but its result could be instable in 
the case the k-th. iteration gives a value of c^ := Ax^ close to zero so that 
the punctual division c q := c./c* is unaccurate. 

We observe from the results that the lowest error for MLEM is err = 0.4260 
at iteration 700, which is obtained after a time t = 4.177 while the lowest 
error for LSCG is err = 0.5668 at iteration 800, which is obtained after a 
time of t = 4.162. 

Now we test our algorithms in the case the Radon data is subject to a 
Gaussian relative noise of null mean and standard deviation equal to a. We 
have chosen the value of a = 10% for the first experiment. As we can see 
in Figure 8.4, that the MLEM algorithm is more accurate but unstable, and 
the LSCG error reflects the typical behavior of the iterative methods in noisy 
cases. In fact the error does not decrease by raising the number of iterations, 
but exists an optimal iteration number, in our case around 200, where the 
error is minimal. This is reasonable, because the iterative algorithm tries to 
reach the image whose sinogram is closer to our noisy data, and it will reach 
a noisy image. The discussion in medical image reconstruction about the 
choice of the optimal iteration remains open, because it is dependent from 
the image itself. 
As a second experiment, we tried to vary the relative noise level from to 

1 it typically shows a good results only after N 2 iterations, i.e. about 16000 for the 
standard case. 

2 http://www2.imm.dtu.dk/^pch/AIRtools/. 
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100%, and use the same number of iteration for both algorithms. In Figure 
8.5 we can observe that the MLEM is more accurate, but again, is "allergic" 
to some noise levels at some iteration and the error of LSCG does not raise 
so much with noise and can bring good results even if the noise level is high. 



" i ii 'i i i I., i . I ? i 



Error for i. iterations of MLEM 
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Time of computation for k iterations of LSCG 



Error for K iterations of LSCG 
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»-norm 

error 



Number of iterations 



Figure 8.3: Time and error test for MLEM (upper) and LSCG (lower) at several 
iterations. 
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Figure 8.5: Error of the MLEM and LSCG algorithms with several noise levels. 
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Figure 8.4: Error of the MLEM and LSCG algorithms with a noise of a = 10% 
after several numbers of iterations. 
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8.3 Kernel methods results 

We made 3 some tests also on the kernel method described in Chapter 5, using 
the four functions we have selected before: the Gaussian function the Ball 
function, the Wendland and the Wu function. As we can see in Figure 8.6, 
the sparsity of the system matrix depends on the shape parameter e: as e 
decreases, the support becomes larger and the system matrix becomes less 
sparse (i.e. it has less elements equal to zero). If the matrix has few nonzero 
elements the results are inaccurate, while on the other side if the nonzero 
elements are too much the time needed for assembling the system matrix 
and the time of solution of the system is very long. Therefore a good trade- 
off between time computation and accuracy is needed. For this reasons we 
have chosen to work on values of e G [0.5, 10]. We made also further studies 
on a cluster with the aim to investigate the results for smallest values of e. 
We have chosen for this test to use the LSCG algorithm for 50 iterations 
because we considered this choice as a good trade-off between speed and 
accuracy. 

We also have found the "optimal" shape parameter with trial-and-error, i.e. 
we have selected the value of e which gives a lower error. Finally we have 
verified the error bound proved in Theorem (5.3.3), i.e. 



11/ - /||l~ < 2\S 1 \ \\w ideal \\ L i ChxflWKfy* 



:(") 



where ||u> i(fea ;|| L i < 6y fg. Now IS" 1 ] = vr, b = ^ = 0.053, then \\w idea i\\ L i < 

^\/T8 ~ O-l^l and hx,n = ^ because the data centers are equispaced with 

unitary distance. 

We have estimated the constant C experimentally, as 

r 11/ "/loo 

C = max 



2\&\bJ ^hxAKfU* 



and verified the error bound. We can see the results of these tests in Figures 
8.7 and 8.8; the constant C the optimal shape parameter e opt are in Table 
8.1. 



3 our scripts are based on the Air Tools package, as all the iterative algorithms in this 
work. 
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Figure 8.6: Sparsity of the classical ART matrix and of the matrix obtained using 
Gaussian kernel and e = 1; the shape parameter e can change the (essential) 
support of the kernel and the sparsity too. A too sparse matrix can be too ill- 
conditioned, while a "richer" matrix slows down the reconstruction. A good trade- 
off is required in the choice of the shape parameter. 



function 


C 


c-opt 


error for e opt 


Gaussian 


0.0406 


1.5 


0.4903 


Ball 


0.0437 


1 


0.4957 


Wendland 


0.0446 


1.5 


0.6011 


Wu 


0.0434 


1 


0.5104 



Table 8.1: This table shows the results of the experiments made on kernel func- 
tions for e G [0.5,10]; for each function is reported the corresponding estimate 
for the constant C, the optimal shape parameter, and the oo-norm error of the 
corresponding reconstruction. 

We made two more tests for e G [0.01, 1]. As we mentioned before we have 
used cluster computing for such values of the shape parameter. To be more 
precise we used the hpbladel4 knot of the NumLab at the Department of 
Mathematics 4 . This knot uses two processors Quad-Core AMD Opteron™, 
2.3GHz, and 64Gb of RAM memory, of which we have used 32Gb. 
In Figure 8.9 we see the results of two different sampling of e: the first one 
is a quasi-equispaced sampling e G {0.01, 0.05, 0.1, 0.15, ..., 1}, the second 
one considers different rays of the (essential) support - G {1, 2, • • • , 20}. As 
we can see from both the results, the optimal shape parameter for Ball and 
Wu functions is still e opt = 1. On the contrary the optimal parameter e opt is 
smaller for the Gaussian (e opt = 0.55) and Wendland (e opt = 0.091). 
Concerning the speed, the iterative methods have the same order of conver- 
gence than the classic iterative methods, O(kN), but are often slower be- 



4 see the page http://numlab.matfi.imipd.it/?q=node/23 for further details. 
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cause they require also the multiplying the coefficients for the data matrix 5 
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Figure 8.7: Time of computation in seconds (left), error and error bound (right) in 
oo-norm of the kernel method for the functions Gaussian (upper) and "ball" (lower) 
with several shape parameters e G [0.5, 10] 



5 if the shape parameter is big it happens that H — I the identity matrix, and the 
computational time required for solving our system is less than the time of computation 
for a classical MLEM or LSCG algorithm. 
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Figure 8.8: Time of computation in seconds (left), error and error bound (right) in 
oo-norm of the kernel method for the functions Wendland (upper) and Wu (lower) 
with several shape parameters e £ [0.5, 10]. 
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Figure 8.9: Results of two tests made with cluster computer for two different 
samplings of e G [0.01, l].The optimal shape parameter e op t is signed with a circle. 



8.4 Comparisons of the results 

We will now compare the results of the different algorithms in terms, as usual, 
of oo-norm error. 

Let us compare the efficiency of the FBP algorithm 6 with the iterative LSCG 
algorithm at 50 iterations and the kernel method using Gaussian function, 
optimal shape parameter and 50 iterations of LSCG method. As we can see 
in Figure 8.10, both the iterative tests lead to more accurate results than 
those of the FBP, and the Gaussian can improve the results of the LSCG. 
If we look at the mean time of computation we see that the FBP runs on 
t\ = 0.1s, while the LSCG and Gaussian methods require to be ran at time 



6 when this time we have forced the dimensions so that it reconstruct a 128x128 image 
from a 128x60 data matrix. 
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t 2 = 0.3s and £3 = 0.8s respectively. This means that the kernel method with 
Gaussian reduces the error of 15% requiring 160% more time with respect to 
the LSCG method. 

As a second comparative test we will vary the number of iterations and 
compare the classic iterative methods with the respective kernel versions, i.e. 
the kernel method using such methods. The comparison will be made on time 
of computation and error. We have chosen for this test the Gaussian function 
which appears to be the most precise and slow of the kernels and for the 
Wendland function, which seems to be the fastest and least accurate. As we 
see in Figure 8.11, the results obtained are what we expected in LSCG case. 
In fact the Gaussian is the slowest but most accurate, while the Wendland 
is the fastest but least accurate, while the classic method offers mean results 
both in terms of speed and of accuracy. To the contrary the MLEM test 
offers a surprise: the Gaussian function's results appear to be faster after a 
few iterations, but the classical method passes the kernel method after 70 
iterations. The Wendland function's results remain always the faster but the 
least accurate between the three methods. 



I 1 1 1 1 1 J 1 1 J ; 1 1 , 777: 



-FBP error 
-LSCG error 
"Gaussian error 
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Figure 8.10: Error of FBP, LSCG, and Gaussian algorithms as the relative error 
varies from to 100%. The results of the FBP formula are worse than those in 
the previous section because we're forcing the FBP algorithm to work on the same 
dimensions of the other algorithms and not on its standard dimensions. 
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Figure 8.11: Computational times in seconds (left) and errors (right) for the LSCG 
and its kernel versions (upper), for MLEM and its kernel versions (lower) with 
Gaussian and Wendland functions. Along x we have the number of iterations. 



73 



8.5 Hybrid SPECT/CT simulations 

In this last section we will compare the results of the two methods we have 
proposed for the hybrid SPECT/CT reconstruction. The first one is the an- 
alytic algorithm based on Novikov-Natter formula 7 , very similar to the FBP, 
while the second one is based on the iterative methods 8 , introduced in Chap- 
ter 6. Since the analytical algorithm does not allow any change of dimension 
we will work on its standards: we will reconstruct a 128x128 image using 
a a data matrix 182x120 and a (reconstructed) attenuation map 128x128. 
For this purpose we are using two phantoms as in Figures 6.2 and 6.3 in 
Chapter 6. The first one is the attenuation map, that will be projected and 
back-projected using the iterative LSCG at 100 iterations. Then we will use 
the reconstructed attenuation map to compute the attenuated Radon trans- 
form of the second phantom and invert the transform using the two methods. 

Before comparing the two methods we will choose the parameter that 
weights the effect of the attenuation. For this purpose we will simply make 
some reconstruction (always using 100 iterations of LSCG) at on some values 
of A G [0, 1] and select the value that gives a smallest oo-norm error. The 
results of this test are shown in Figure 8.12. As we can see the error can 
change very much, while the total time needed to compute the system matrix 
and to solve the system is almost constant. This test indicates us the value 
of X op t = 0.04 as the "optimal" value, which gives an error of e = 0.235. 

Now we can compare the two algorithms in terms of oo-norm error of 
the reconstruction of noisy signals. As we have done in previous section we 
will compare the errors of the two methods when the data are affected from 
a relative noise which goes from to 100%. For the iterative method we 
have used again 100 iterations of LSCG algorithm and the optimal value 
^opt = 0.04 found in the previous experiment. 

From Figure 8.13 we can see that the time required for the solution of the 
system with iterative method is a little longer than the time used by the 
analytical one. On the other hand the iterative formula reveals a lowest 
error and a better resistance to the noise. 



7 for our tests we have used the functions made by F. Monard and available at his 
webpage http://www.math.washington.edu/~fnionard/. 

8 and then is based on the Air Tools package as all the iterative algorithms in this work. 
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Figure 8.12: Trial-and-error test for the A parameter. To the left we see the total 
computational time of the matrix B and the solution of the system. Right: the 
errors for various values of A £ [0, 1]. 
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Figure 8.13: Time of resolution in seconds (left) and error (right) for the analyt- 
ical and iterative methods for the resolution of the hybrid SPECT/CT simulated 
problem at several relative noise levels form to 100%. 
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Chapter 9 

Results from In vitro experiments 



As far as the laws of mathematics refer to reality, they 
are not certain; and as far as they are certain, they do 
not refer to reality. 

Albert Einstein, quote reported in the 1960 
article The Unreasonable Effectiveness of Mathematics in 
the Natural Sciences by E. Winger. 

)p"?W™^(\ his chapter is dedicated to some experiments made on duin- 
S^?!/^ niies (in vitro). The quality of the results is not evaluated 
^^Bf^S' scientifically, we will only make some visual considerations. 
^zM!!^^ We will also discuss about the resolution of the Gamma camera 
we are using, and compare it with the experimental resolution. 

9.1 Resolution and sampling 

Before we discuss about the laboratory results, we will verify that the stan- 
dard sampling (120x128) is correct with respect to the system resolution and 
does not give data with lower resolution than the resolution of the machine. 
First of all we have to compute the machine resolution, i.e. the resolution of 
the collimator in use in our machine. 

The gamma camera used in the experiments is the Ultra High Resolution 
(UHR), and the collimators have a hexagonal structure, like the hives of 
bees, and the crystal type is 3/4". 

In Table 9.1 are reported the main characteristics of the most commons ma- 
chines used in SPECT. 
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Description 


Max. 
useful 

energy 
(Kev) 


Septal 
thick- 
ness 
(mm) 


Hole size 
(mm) 


Hole 

length 

(mm) 


Resolution 
at 0mm 
(FWHM, 
mm) 


Resolution 
at 100mm 
(FWHM, 
mm) 


Crystal type 










3/8"; 3/4" 


3/8"; 3/4" 


Dynamic 


150 


0.305 


2.54 


25.4 


5.1; 5.7 


14.8; 15.0 


High Resolu- 
tion 


160 


0.203 


1.22 


27.0 


3.9; 4.7 


7.0; 7.6 


UHR 


160 


0.152 


1.78 


58.4 


4.1; 4.8 


6.0; 6.7 


Medium en- 
ergy 


300 


0.864 


3.40 


58.4 


4.1; 4.8 


6.0; 6.7 


High energy 


440 


2.01 


3.4 


58.4 


5.9; 6.6 


11.0; 11.4 



Table 9.1: Main characteristics of the most commons machines used in SPECT. 



The declared intrinsic resolution is Ri = 4.1mm. 
Let us take a look on the septal thickness. From our formula we find 



t> 



6D 



6 • 1.78mm 



jj,l — 3 2.14mm _1 • 58.4mm — 3 



0.088mm 



The septal thickness of our machine, t = 0.152mm, fulfills the condition 

above. 

The collimator resolution is R = D(l + t^—) if we consider that the collima- 

tors are made of lead, \x = 2.14mm _1 at 150 keV, then the effective collimator 

length is 

2 

L e ff = L = 58.4mm — 0.9mm = 57.5mm 

/i 

So we can compute the collimator resolution 



R r 



1.78 
57^5 



x + 1.78 = 0.031a; + 1.76 



This means that the resolution at mm is 1.78mm and at 100 mm we have 

a resolution of 4.9 mm. 

We can also calculate the sensitivity 

, , KD\ 2 /0.26 • 1.78mm\ 2 „ ,„ ,- 
S=[ = ^6.47-10 -5 . 



J eff 



57.5mm 
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What is the minimal sampling which guarantees a resolution of R s = 6.7mm 
(the biggest resolution at 100mm)? Now, this resolution means a maximum 

frequency at 

2vr 
b = —100 =* 93.78 
R s 

so the minimum number of angular samples is / = \b\ + 2 = 95 and the 
minimum number of equispaced linear samples is p = |_— J = 59. This means 
that the standard sampling is correct and does not decrease the resolution. 
We can also obtain from this data an estimate of the system resolution, which 
is R s = s/Rf + R%, where we know the value of Ri = 4.1mm. This means 
that the system resolution is R s = 4.47mm at 0mm and R s = 6.8mm at 
100m, which is consistent with the values of FWHM declared on Table 9.1. 



9.2 Three capillaries experiment and experi- 
mental resolution 




Figure 9.1: The three capillary dummy used in experiments. 

The first dummy used in experiments, shown in Figure 9.1, is a cylinder 
filled with water with three thin tubes on the inside, filled with tracer. This 
experiment is a valid alternative to the small sphere filled with tracer for 
the estimation of the experimental resolution as FWHM of the reconstructed 
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image. The data are taken with different radius of acquiring, precisely r\ = 

145mm, r 2 = 165mm and r 3 = 185mm. With this data we will be able to 

extrapolate an estimation of the resolution at a desired distance. 

The raw data is a 128x128x120 matrix which contains the projection data 

taken in 120 angular projection by a grid of 128x128 collimators. This is 

equivalent for our purposes to 128 sinograms 128x120, each one of them 

represents the analysis at different heights. 

In Figure 9.2 we can see a reconstruction and segmentation of a slice of this 

dummy, obtained using the LSCG at 100 iterations and the phantomization 

algorithm. 




iterative reconstruction at height 86 on 123. 



segmentation of the reconstruction 





Figure 9.2: Raw data sinogram, reconstruction and phantomization for the three 
capillaries dummy data. 

From the reconstruction of the central capillary at different slices we 
obtain some estimates of the FWHM 1 . The mean values of those FWHM 
are reported on Table 9.2. 



distance 


145mm 


165mm 


185mm 


FWHM 


8.78 


9.53 


10.62 



Table 9.2: Values of the experimental FWHM at several distances between abject 
and collimators. 

From these three values we can extrapolate which could be a reasonable 
experimental resolution at a given distance. In fact the intrinsic resolution 
Ri = 7 is a constant, while the collimator resolution depends linearly from 
the distance R c = ax + (3, with a, 7, fi > 0. So the system resolution will be 



1 this work has done by prof. Cecchin, we want to thank him again for his fruitful 
cooperation. 
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of the form 



FWHM{x) =R S = ^R 2 + R 2 = ^a 2 x 2 + 2a[5x + (3 2 + 7 2 = Vax 2 + bx 

with a,b,c > 0. 

If we want to find the constants a,b,c we have to solve this problem 

min || a/ox 2 + bx + c — f\\ 2 

a,b,c 



with the constraints a,b,c > 0. In our notation now x = (45, 165, 185) the 
distance vector and / = (8.78, 9.53, 10.62) are the data, and the chosen norm 
is the discrete 2-norm. 

To solve this problem we have used the function fmincon contained in the 
MATLAB® package Optimization Toolbox. The results are shown in Figure 
9.3. From these results we find a resolution of 4.71mm at 0mm and of 6.79mm 
at 100mm, which is close to the declared values of FWHM in Table 9.1. 



experimental fwhm 
" extrapolation 



20 40 



Figure 9.3: Extrapolation results from FWHM data using a nonlinear optimization 
function. 
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9.3 Cerebral dummy 




Figure 9.4: The cerebral dummy used in experiments. 

The second dummy used in experiments is a cerebral dummy like the one 
showed in Figure 9.4. This dummy is done with materials with the same 
attenuation coefficients of the human head. Instead of the brain is filled with 
water (the attenuation coefficients are similar) and tracer, while the four 
buffers on the inside are filled with a more concentrated tracer. 
Again the raw data is a 128x128x120 matrix with 128 slices at different 
heights. In this case the data is very "poor", because the matrix contains only 
integral numbers between and 17, while the previous dummy data were up 
to 39. The low variance of the data is an index of a great attenuation and 
low concentration of the tracer 2 . So it is more difficult to reconstruct this 
data than the numerical phantom. As said before, we will only show some 
reconstructions and comment them. 



2 the tracer concentration and activity are the same used in real analysis on patients, 
so it's very low for security reasons. 
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sinogram at height=92 on 1 28. 



analytical reconstruction 





iterative reconstruction 



Kernel reconstruction 





Figure 9.5: Different reconstructions of the sinogram data at slice number 92. 

We have selected one slice and we have made different reconstructions, 
shown in Figure 9.5, and we have put all the data less than a chosen 5 to 
zero to eliminate some noise. The first reconstruction is made with FBP, 
5 = 4. We have also added some zeros to the data for reaching the standards 
dimensions of the algorithm As we see from the Figure, we can recognize the 
buffers, but the shape of the dummy is not clear. The second reconstruction 
method is an iterative MLEM after 100 iterations. We can identify both the 
buffers and the dummy shape in this image. The last method is a kernel 
interpolation using the Gaussian function with e — 1.5 and an MLEM at 100 
iterations. In this case we can recognize clearly the shape of the dummy, but 
the buffers on the inside are disappeared. 
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Conclusions and future work 



Briefly, by points, we will summarize the work done and what remains to do. 

Work done 

• Introduced CT and SPECT tomographies both from a medical and 
from a physical point of view (see Ch. 1). 

• Studied the state-of-art algorithms for CT and SPECT (see Ch. 3-4). 

• Found an error bound for the FBP formula (see Ch. 3.4). 

• Introduced Kernel Methods with a relative error bound (see Ch. 5). 

• Proposed an iterative method for SPECT/CT (See Ch. 6.4). 

• Studied some methods for the regularization (see Ch. 7). 

• Made some In silico and in vitro experiments (See Ch 8-9). 

• Studied the connection between the resolution of the machine, of the 
sampling and experimental (see Ch 1.6, 3.2, 9.1). 

• Proposed a method of extrapolation of the experimental resolution at 
a given distance from some data (see Ch 9.2). 

To do: 

• Find an error bound for SPECT/CT algorithms 

• Implement OSEM algorithm and test it with Kernel methods 

• Study the setting of A in iterative SPECT/CT 

• Find an algorithm for contour detection in "poor-data" case 
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Further in vitro experiments for the Kernel methods 
Extend the Kernel methods for 3D reconstructions 
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